Exercise 5-2: Oscillations in a Vertical Plane on Three PlatformsΒΆ
Here particles oscillate through the center within a vertical plane. The period for each platform is displayed once it finishes.
Index 0 is simple harmonic motion - Index 1 is a paraboloid - Index 2 is a concave-up hemisphere
Code translated from GW-BASIC provided in Exercise 5.2 of Stommel and Moore (1989)
"""
author: Victoria McDonald
email: vmcd@atmos.washington.edu
website: https://github.com/torimcd/coriolis-sm
"""
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from matplotlib.animation import FuncAnimation
import numpy as np
from IPython.display import HTML
%matplotlib inline
Here we set up the constants and initialize the global variables
# set up global variables
dt = 0.01 # time step
t = 0
g = 1 # gravity
d = 1
k = 1 # spring constant
k2 = k**2
c = k2/2
a = k2/g
mu = 0 # angular momentum
r = {}
dr = {}
ddr = {}
p = {}
dp = {}
x = {}
y = {}
flag = {0:0, 1:0, 2:0}
drl = {0:-1, 1:-1, 2:-1}
zzz = {}
x_plane = []
y_plane = []
x_para = []
y_para = []
x_hemi = []
y_hemi = []
l_plane = []
l_para = []
l_hemi = []
m_plane = []
m_para = []
m_hemi = []
# set up the surfaces
for i in range(3):
r[i] = 0.5
dr[i] = 0
p[i] = 0
Here we set up the figure that the particles will travel on
%%capture
fig, axes = plt.subplots(2,3)
# turn the axes off
axes[0,0].axis('off')
axes[0,1].axis('off')
axes[0,2].axis('off')
axes[1,0].axis('off')
axes[1,1].axis('off')
axes[1,2].axis('off')
# set the aspect ratio
axes[0,0].set_aspect('equal')
axes[0,1].set_aspect('equal')
axes[0,2].set_aspect('equal')
# create outline
outline_1 = patches.Ellipse((0,0), 1, 1, 0, linewidth=2, fill=False)
outline_2 = patches.Ellipse((0,0), 1, 1, 0, linewidth=2, fill=False)
outline_3 = patches.Ellipse((0,0), 1, 1, 0, linewidth=2, fill=False)
# add outline to axes
axes[0,0].add_patch(outline_1)
axes[0,1].add_patch(outline_2)
axes[0,2].add_patch(outline_3)
axes[0,0].set_xlim(-0.6, 0.6)
axes[0,0].set_ylim(-0.6, 0.6)
axes[0,1].set_xlim(-0.6, 0.6)
axes[0,1].set_ylim(-0.6, 0.6)
axes[0,2].set_xlim(-0.6, 0.6)
axes[0,2].set_ylim(-0.6, 0.6)
axes[1,0].set_xlim(-1, 1)
axes[1,0].set_ylim(-1, 1)
axes[1,1].set_xlim(-1, 1)
axes[1,1].set_ylim(-1, 1)
axes[1,2].set_xlim(-1, 1)
axes[1,2].set_ylim(-1, 1)
plt.suptitle("Three Surfaces of Revolution")
axes[1,0].text(-0.5, -0.5, "plane")
axes[1,1].text(-0.5, -0.5, "paraboloid")
axes[1,2].text(-0.5, -0.5, "hemisphere")
Now we set up the particle objects to be updated in each step of the animation
# set up the particles
plane = axes[0,0].plot([], [], color='green')[0]
paraboloid = axes[0,1].plot([], [], color='red')[0]
hemisphere = axes[0,2].plot([], [], color='gold')[0]
plane_section = axes[1,0].plot([], [], color='green')[0]
paraboloid_section = axes[1,1].plot([], [], color='red')[0]
hemisphere_section = axes[1,2].plot([], [], color='gold')[0]
This is the update() function that calculates the new position for each surface, and updates the particle objects at each step
n=500
def update(frame):
global dt
global g
global d
global k
global k2
global c
global a
global mu
global r
global dr
global ddr
global p
global dp
global x
global y
global t
global x_plane
global y_plane
global x_para
global y_para
global x_hemi
global y_hemi
global l_plane
global l_para
global l_hemi
global m_plane
global m_para
global m_hemi
global flag
global drl
# draw cross sections first
if frame < 70:
i_map = np.arange(-0.7, 0.7, 0.02)
i = i_map[frame]
l_plane.append(i)
m_plane.append(0)
l_para.append(i)
m_para.append(c*(i**2))
if (a**2 - i**2) <= 0:
plane_section.set_data(l_plane[:], m_plane[:])
paraboloid_section.set_data(l_para[:], m_para[:])
hemisphere_section.set_data(l_hemi[:], m_hemi[:])
else:
l_hemi.append(i)
m_hemi.append(a-np.sqrt(a**2 - i**2))
plane_section.set_data(l_plane[:], m_plane[:])
paraboloid_section.set_data(l_para[:], m_para[:])
hemisphere_section.set_data(l_hemi[:], m_hemi[:])
else:
# integrations for plan view
ddr[0] = mu**2/(r[0]**3)-k**2*r[0]
aaa = (1+(2*c*r[1])**2)
ccc = mu**2/(r[1]**3)
bbb = -2*g*c*r[1] - r[1]*(2*c*dr[1])**2
ddr[1] = (bbb+ccc)/aaa
aa = -g*r[2]/(np.sqrt(a**2 - r[2]**2))
bb = mu**2/(r[2]**3)
cc = -dr[2]**2*(a**2*r[2]/((a**2-r[2]**2)**2))
dd = 1+r[2]**2/(a**2*r[2]**2)
ddr[2] = (aa+bb+cc)/dd
for i in range(3):
dr[i] = ddr[i]*dt + dr[i]
r[i] = dr[i]*dt + r[i]
#for j in range(3):
dp[i] = mu/(r[i]**2)
p[i] = dp[i]*dt + p[i]
#for k in range(3):
x[i] = r[i]*np.cos(p[i])
y[i] = r[i]*np.sin(p[i])
#for l in range(3):
if flag[i] == 1:
drl[i] = dr[i]
else:
zzz = dr[i]*drl[i]
if zzz<0:
flag[i] = 1
axes[1, i].set_title("time("+str(i+1)+"): " + str(2*t)[0:5], loc='left')
else:
drl[i] = dr[i]
t = t+dt
x_plane.append(r[0]*np.cos(p[0]))
y_plane.append(r[0]*np.sin(p[0]))
x_para.append(r[1]*np.cos(p[1]))
y_para.append(r[1]*np.sin(p[1]))
x_hemi.append(r[2]*np.cos(p[2]))
y_hemi.append(r[2]*np.sin(p[2]))
plane.set_data(x_plane[:], y_plane[:])
paraboloid.set_data(x_para[:], y_para[:])
hemisphere.set_data(x_hemi[:], y_hemi[:])
Finally, we construct the animation and embed it in this page
# Construct the animation, using the update function as the animation director.
animation = FuncAnimation(fig, update, frames=n, interval=10, repeat=False, blit=False)
# convert to a video to be embedded in web page
HTML(animation.to_jshtml())